---
title: "ASE"
output: html_document
---

```{r setup, include=FALSE}
library(DESeq2)
library("magrittr")
library("dplyr")
library("ggplot2")
library("reshape2")
library("gplots")
library("ASEP")
```

#Read in data
```{r cars}

phased = read.table("Edited_read_counts_for_ASE.csv", sep=",", header=TRUE)

```

#single condition
```{r pressure, echo=FALSE}


genes = unique(phased$gene)


# Run ASEP, note that because of resampling P-values will be slightly different every time

ase_table=rep(NA,1)

for (i in seq(1,length(genes))) {
  
gene1 = genes[i] 

data = subset(phased, gene==gene1)

new = ASE_detection(dat_all=data, phased=TRUE, adaptive=TRUE, n_resample=10^3, parallel=FALSE, save_out=FALSE)

print(new)
rbind(ase_table, new)

}


#Read in analysis used in the paper and do an FDR
ase.out = read.table(file="ASE_out.csv", sep=",", header=TRUE)
ase.out$Pvalue = as.numeric(ase.out$Pvalue)
ase.out$fdr = p.adjust(ase.out$Pvalue, method = "fdr", n = length(ase.out$Pvalue))


ase_sig = subset(ase.out, fdr < 0.1) #FDR cutoff of 0.1

sig = unique(ase_sig$gene) #Get list of significant genes


```

#Run topGO
```{r}
library(topGO)
library(splitstackshape)

#read in list of contigs and whether they show alpha biased expression, beta biased expression, or allele biased expression
to_test = read.table("for_GO.csv", sep=",", header = TRUE)

#read in GO terms associated with each gene
geneID2GO = readMappings("ase_go.txt")


#Run topGO
geneNames = names(geneID2GO)

#subset out the three groups
ase_a = to_test$transcript_id[to_test$a_exp==1]
ase_b = to_test$transcript_id[to_test$b_exp==1]
ase_u = to_test$transcript_id[to_test$ase_exp==1]


#alpha biased expression

geneList = factor(as.integer(geneNames %in% ase_a))
names(geneList) <- geneNames



GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,
                    annot = annFUN.gene2GO, gene2GO = geneID2GO)

resultTopGO.elim <- runTest(GOdata, algorithm = "elim", statistic = "Fisher" )
allRes_alpha <- GenTable(GOdata, elimKS = resultTopGO.elim,
                   orderBy = "elimKS", 
                   topNodes = length(geneList))

allRes_alpha$padjust<-p.adjust(allRes_alpha$elimKS,method="BH")


#beta biased expression

geneList = factor(as.integer(geneNames %in% ase_b))
names(geneList) <- geneNames



GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,
                    annot = annFUN.gene2GO, gene2GO = geneID2GO)

resultTopGO.elim <- runTest(GOdata, algorithm = "elim", statistic = "Fisher" )
allRes_beta <- GenTable(GOdata, elimKS = resultTopGO.elim,
                   orderBy = "elimKS", 
                   topNodes = length(geneList))

allRes_beta$padjust<-p.adjust(allRes_beta$elimKS,method="BH")

#Allele biased expression 

geneList = factor(as.integer(geneNames %in% ase_u))
names(geneList) <- geneNames



GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList,
                    annot = annFUN.gene2GO, gene2GO = geneID2GO)

resultTopGO.elim <- runTest(GOdata, algorithm = "elim", statistic = "Fisher" )
allRes_u <- GenTable(GOdata, elimKS = resultTopGO.elim,
                   orderBy = "elimKS", 
                   topNodes = length(geneList))

allRes_u$padjust<-p.adjust(allRes_u$elimKS,method="BH")

```



#plot Figure 4 and Supp Fig 12
```{r}

#Read in data for Supp Fig 12

all = read.table("all_data_forgraph.csv", sep=",", header=TRUE) #read in data for all loci
sig = as.data.frame(sig)
sig$sig = gsub("_TRINITY_", "_", sig$sig) #Change names of genes to match so that they fit in the headers
all$sig = match(all$gene, sig$sig, nomatch = 0)
all$sig[all$sig > 0] <- 1 

all$id = as.factor(all$id)
all$sig = as.factor(all$sig)

all_sig = subset(all, sig==1)

#Read in data for Figure 4

new = read.table(file="sig_data_graph.csv", sep=",", header=TRUE)

new$increased_exp = factor(new$increased_exp, c("A", "U", "B"))



#Special function that allows the 1:1 line be the same in all facets

FacetEqualWrap <- ggproto(
  "FacetEqualWrap", FacetWrap,
  
  train_scales = function(self, x_scales, y_scales, layout, data, params) {
    
    # doesn't make sense if there is not an x *and* y scale
    if (is.null(x_scales) || is.null(x_scales)) {
        stop("X and Y scales required for facet_equal_wrap")
    }
    
    # regular training of scales
    ggproto_parent(FacetWrap, self)$train_scales(x_scales, y_scales, layout, data, params)
    
    # switched training of scales (x and y and y on x)
    for (layer_data in data) {
      match_id <- match(layer_data$PANEL, layout$PANEL)
      
      x_vars <- intersect(x_scales[[1]]$aesthetics, names(layer_data))
      y_vars <- intersect(y_scales[[1]]$aesthetics, names(layer_data))
      
      SCALE_X <- layout$SCALE_X[match_id]
      ggplot2:::scale_apply(layer_data, y_vars, "train", SCALE_X, x_scales)
      
      SCALE_Y <- layout$SCALE_Y[match_id]
      ggplot2:::scale_apply(layer_data, x_vars, "train", SCALE_Y, y_scales)
    }
    
  }
)


facet_wrap_equal <- function(...) {
  # take advantage of the sanitizing that happens in facet_wrap
  facet_super <- facet_wrap(...)
  
  ggproto(NULL, FacetEqualWrap,
    shrink = facet_super$shrink,
    params = facet_super$params
  )
}

#Make Supplemental Fig 12 
pdf("Sup_fig_12.pdf",width=12,height=10)
ggplot(all_sig, aes(x=ref,y=alt)) + geom_point(size=2) +  geom_abline(slope=1) +  theme(legend.position = "none", axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size=12), axis.text.y = element_text(size=12), strip.text.x = element_text(size=7)) + labs( y = "Read Counts"~paste(beta), x="Read Counts"~paste(alpha)) + facet_wrap_equal(~gene,scales= "free")
dev.off()


#make Figure 4
pdf("Figure4.pdf",width=10,height=10)
ggplot(new, aes(x=ref,y=alt, color=increased_exp)) + geom_point(size=2) + theme_bw() + geom_abline(slope=1) +  theme(legend.position = "bottom", axis.title.x = element_text(size=14), axis.title.y = element_text(size=14), axis.text.x = element_text(size=10), axis.text.y = element_text(size=10), strip.text.x = element_text(size=8), legend.title = element_blank()) + labs( y = "Number of"~paste(beta)~"reads", x="Number of"~paste(alpha)~"reads") + facet_wrap_equal(~gene,scales= "free") + scale_color_hue(labels = c(paste(alpha)~"biased expression","allele biased expression",paste(beta)~"biased expression"))
dev.off()

```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
